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The stability of dark solitons generated by a supersonic flow of a Bose-Einstein condensate past 
a concave corner (or a wedge) is studied. It is shown that solitons in the dispersive shock wave 
generated at the initial moment of time demonstrate a snake instability during their evolution to 
stationary curved solitons. Time of decay of soliton to vortices agrees very well with analytical 
estimates of the instability growth rate. 
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Flow of quantum liquids has attracted much attention since discovery of superfluidity. At first it was discussed 
mainly in relation with the existence of the critical flow velocity past obstacles (or boundaries of capillary tubes) above 
which superfluidity is lost. According to Landau criterium it occurs at the velocity corresponding to the threshold of 
generation of linear excitations in the fluid (for instance, rotons in Hell). However, it was found that the experimental 
value of the critical velocity in Hell was much less than that predicted by Landau and this contradiction between the 
theory and the experiment was explained by Feynman who indicated that the threshold for generation of nonlinear 
excitations, namely vortices, is much less than the Landau value and agrees, at least semi-quantitatively, with the 
experimental value. Later these predictions were confirmed by analytical and numerical study of a superfluid flow 
| modeled by the Gross-Pitaevskii equation for a weakly interacting Bose gas [H, 0. As was found in these studies, a 
p I . 2D flow past a disk-shaped obstac le with the radius not much greater than the healing length looses its superfluidity 
at velocity equal to v cr = yj2/\\c s = 0.43c s , c s being the sound velocity, i.e. the velocity of Bogoliubov excitations 
in the long wavelength limit. Here v cr equals to the threshold velocity above which vortices are generated by the flow 
at the boundary of the obstacle. It was noticed also in the numerical simulations that a stationary wave pattern is 
generated by a supersonic flow past an obstacle and this pattern consists of two different regions separated by the 
Mach cone. These theoretical findings got new interest after experimental realization of the flow of Bose-Einstein 
condensate (BEC) of dilute atomic gas past obstacles Q. In theoretical works motivated by these experiments it was 
shown that the region outside the Mach cone is occupied by the so called "ship waves" formed due to interference of 
Bogoliubov waves radiated by the obstacle by means of Cherenkov effect > whereas in the region inside the Mach 
cone a fan of oblique dark solitons is located @, 0] and these oblique solitons are effectively stable if the flow velocity 
' is large enough @. 

However, if we consider flow of BEC past an extended obstacle, the wave pattern can change drastically. In 
' particular, the ship wave not far from the obstacle's boundary becomes nonlinear and transforms into a stationary 
dispersive shock wave, i.e., in some approximation, into a lattice of curved dark solitons [![. In vicinity of the boundary 
of the obstacle both the local sound velocity and the direction of the flow are changed compared with their values in 
the incident flow and, hence, the question of stability of dark solitons in the dispersive shocks should be reconsidered. 
Their instability was noticed in numerical simulations of the highly supersonic flow past a slender wedge Q and it 
became crucially important for moderately supersonic flows past not so slender wedges when the instability destroyed 
very fast the dispersive shock pattern [Hj. The aim of this paper is to give theoretical estimates of the instability 
growth rate in such flows and to confirm the theory by numerical simulations. 

First, we notice that if the obstacle is large and not too slender, then transition to a stationary state formed by the 
"curved solitons" is accompanied by propagation of bending waves along such solitons; these waves include unstable 
modes which evolution can lead eventually to decay of oblique solitons to vortex-antivortex pairs. Let at the initial 
moment of time t = the condensate be confined in the region outside the trap potential having a form of a corner 
with one its wall located along negative x axis (— oo < x < 0) and another one along the line y = tana • x, a 
being the angle of the slope (0 < a < tt/2). The flow with velocity M (equal to the Mach number in our standard 
non-dimensional units with the flow velocity measured in units of the sound velocity) along the x axis is switched on 
at t = 0. Evidently, in vicinity of the sloping boundary the wave appears which is equivalent to one generated by 
a piston moving with velocity Msina; hence, a non-stationary dispersive shock wave is generated [9l4ll| with wave 
crests parallel to the sloping boundary. If Msina < -,/no, no being the density of incident BEC, then the waves 
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FIG. 1: Density plot of the wave pattern generated by the flow of BEC with the Mach number M — 2.5 past a corner with the 
angle a — 52° at the moment of evolution t = 8.22. The oblique solitons of the dispersive shock are modulated by the wave of 
bending. These wave packets include the unstable modes which lead eventually to decay of oblique solitons to vortices. 



closest to the boundary can be considered as oblique solitons [7|. However, these oblique solitons are cut at their edge 
located about the line y = — cot a ■ x. These their edges are not attached to the obstacles and, hence, cannot become 
effectively stable due to transition from their absolute instability to a convective one at large enough flow velocity (see 
8[). Therefore these oblique solitons must change their form transforming to "curved solitons" studied analytically in 
9] for M 3> 1, a ^ 1. Hence, bending waves start their propagation along oblique solitons and we suppose that the 
component of the wave packet with maximal value of the instability growth rate T leads to decay of oblique solitons 
to vortices observed numerically in [l(3] ■ 

To perform quantitative comparison of analytical estimates with results of numerical simulations, we define a 
numerical time T(a) of the instability growth as follows. Direct numerical solution of the GP equation 

ii>t + \(^xx + ^yy)- \^\ 2 i> = V( X ,y)ljj (1) 

with the corner-like shape of the trap potential V{x, y) shows that solitons bend in the shock with time in such a way 
that, for instance, the deepest one is centered around the curve y' — Y(x',t), where x' is the axis along the sloping 
side of the corner and y' is the axis normal to x' , that is (x',y') axes are obtained by counter-clockwise rotation of 
(x, y) axes to the angle a and Y(x' , t) has a meaning of the distance of the soliton's trough line from the boundary at 
the moment t. A typical shape of this wave pattern is shown in Fig. 1 and one can see that a wave packet propagates 
along the soliton (x' axis). Our calculations show that the packet's velocity is equal to M cos a, that is this packet 
is built from non-propagating unstable modes and its motion is caused by convection due to flow of BEC along the 
boundary. In Fig. 2 the forms of wave packets corresponding to different values of a at such moments of time that 
their amplitudes reach a certain value chosen the same one for all a. As we see, the forms of these wave packets are 
very similar to each other, too, and we choose just these moments of time as definitions of numerical values of the 
instability growth time T(cv). 

To estimate the analytical growth rate values, we notice that the problem of formation of dispersive shock wave by 
a piston moving with a constant velocity is equivalent to well-known problem of evolution of a step-like pulse [H, 03 
and its solution was discussed in 0, [Tl| - In particular, it was found that the condensate's density between the corner's 
boundary (piston) and the soliton edge of the dispersive shock wave is equal to 

n = (y/na + ^Msina) 2 , (2) 

where no is the density of BEC in the incident flow. The velocity of the deepest soliton at the edge of the shock is 
equal to 

V = y/no~+ |Msina. (3) 

These formulae determine the parameters of the soliton which stability must be studied. Its analytical form is well 
known and can be expressed as 

V^ = $ s (C)e- mt , *.(0=Atanh(AC)-^=VnI(Oe i *" {C) , ^ = -V ( 1 - - 



(4) 



\Jn- V 2 , ( = x + Vt. 



FIG. 2: Line of trough of the deepest dark soliton generated by the flow of BEC with M — 2.5 at the moment of time t = 8.22 
past a corner with the angle a — 52°. In the insert the wave packets with the same amplitude are shown for the angles (1) 
a = 52°, (2) a — 37°, (3) a = 32°. The corresponding moments of time are chosen as times T(a) of the instability growth. 



Then the perturbed wave function reads 



t/j = ($ a + cj>)e- int = ($ s + <f>' + i<p")e- int = ^n s + 5ne l( - v ° +Slp) - mt , (5) 
where the real variables <fJ and <j>" are related with perturbations of the density Sn and the phase 5(f) by the formulae 

, cosip s . „ smtp s — 

on + yris sm ip s ■ dip, <p — on — ^Jn s cos ip s - dip. (6) 



2^n~ s v r. r , t 2 

Substitution of ([5]) into ([1]) and linearization of the resulting equation with respect to <f> gives 

ifa +n( j )+ i (4, xx + <p yy ) - (2|$ s | 2 + $^*) = o, (7) 

where asterisk denotes complex conjugation. We are interested in propagation of a harmonic disturbance of the 
soliton along y axis. As usual in linear problems, such a disturbance can be represented in a complex form, 8n, dip cx 
exp(Ti + ipy) or <j) = (<j)' + i(f>") cxp(rt + ipy), <j)* = {<p' — icf)") exp(ri + ipy), where we keep previous notation for the 
amplitudes <f>' — <f>'(Q, 4>" — 4>"(Q- For this form of the disturbance, Eq. ([7]) reduces to 

iTcf, + iVcf> c + n<p + i0 cc - \y 2 <\> - (2|$ s | 2 <?!> + ^\<f) = 0. (8) 

Substitution of (U]) and separation of real and imaginary parts yields the system 

where 

L n = -Vd c - 2fcT/tanh(fcC), £12 = ~kd CC -n + n s + 2V 2 + \p 2 , 

(10) 

L 2 i = ±d c( + n - 3n s + 2V 2 - \p 2 , L 22 = -Vd ( + 2fcFtanh(fcC). 

This eigenvalue problem was studied in [14] and we used it for numerical calculation of the instability growth rate 
F(p) with the background density n and the soliton's velocity V given by ((2|) and ([3]), respectively, for M — 2.5 and 
several values of the angle a. These values of a as well as corresponding values of the background density and solitons 
velocities are given in the first three columns of Table 1. The instability growth rate depends on the wave number 
p and its maximal value T max for each a is indicated in the fourth column of Table 1. As we see, it decreases with 
growth of a; for convenience we have normalized its values dividing them by the value at a = 52°: 

/ \ r maa; (a) ftt\ 

JtheorW = = tt^t- (11) 

t max \OZ J 

In the last two columns we compare these theoretical predictions with the corresponding numerical characteristics 
defined as 

T(52°) 

Inum(a) = -^—^ (12) 

T[a) 



where T(a) is the defined above instability growth time. One can see very good agreement which confirms our 
approach to the problem of instability of oblique solitons. 
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Table I. Analytical and numerical estimates of instability growth rate of oblique solitons (Mach number M = 2.5). 
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Thus, our results show that the size of the obstacle is crucially important for realization of transition to effectively 
stable dark soliton states. In numerical experiment of Ref. 0] and theoretical estimates of Ref. || the size of the 
obstacle did not play any essential role. However, if a characteristic size I of the obstacle is such that the time l/v conv 
(where v conv is the velocity of convection) is greater than the inverse instability growth rate , then the unstable 
wave packet has enough time for evolution into its nonlinear stage with formation of vortex-antivortex pairs. Hence, 
if we wish to observe effectively stable dark solitons in the flow of BEC past a concave corner or a wedge, then the 
parameters of the flow must satisfy the condition 

; « ^— (13) 



This condition was fulfilled in most of numerical experiments of Ref. |9j but was not fulfilled in numerical experiments 
of Ref. [13] what has lead to fast decay of oblique solitons into vortices. Our estimates show that the location of the 
"seed" disturbance where the first vortices are formed is determined by the distance M cosa/T max along the inclined 
wall of the corner. 

It is worth noticing also that importance of existence of such a period in evolution of dark solitons that the unstable 
modes have enough time for their development into vortices manifested itself also in evolution of dark ring solitons 
[TBI [l6j : the instability develops effectively near the turning point of the ring soliton evolution when its parameters 
do not practically change and the mode with maximal corresponding value of T dominates in the evolution of the 
soliton's shape. x 
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